arXiv: 1501.0209lvl [cond-mat.quant-gas] 9 Jan 2015 


Dipolar Bilayer with Antiparallel Polarization — a Self-Bound 

Liquid 

Martin Hebenstreit 1,2 , Michael Rader 1,2 , and R„ E. Zillich 2 
1 Institute for Theoretical Physics, University of Innsbruck, 

Technikerstr. 25, 6020 Innsbruck, Austria 
2 Institute for Theoretical Physics, Johannes Kepler University, 
Altenbergerstrasse 69, fOfO Linz, Austria 

Abstract 

Dipolar bilayers with antiparallel polarization, i.e. opposite polarization in the two layers, exhibit 
liquid-like rather than gas-like behavior. In particular, even without external pressure a self-bound 
liquid puddle of constant density will form. We investigate the symmetric case of two identical 
layers, corresponding to a two-component Bose system with equal partial densities. The zero- 
temperature equation of state E(p)/N, where p is the total density, has a minimum, with an 
equilibrium density that decreases with increasing distance between the layers. The attraction 
necessary for a self-bound liquid comes from the inter-layer dipole-dipole interaction that leads 
to a mediated intra-layer attraction. We investigate the regime of negative pressure towards the 
spinodal instability, where the bilayer is unstable against infinitesimal fluctuations of the total 
density, conformed by calculations of the speed of sound of total density fluctuations. 
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Introduction. Experiments with Bose gases of atoms with large magnetic moments 
( 52 Ci®21, 164 Dy^, 168 Er,S) are fueling the interest to understand the effects of the dipole- 
dipole interaction (DDI) on stability, shape, and dynamics of dipolar Bose condensates 
(reviewed in Refs. 5 '). The strength of the DDI can be characterized by the dipole length 
r'o = mD 2 /(iiTEoh 2 ), where m is the mass of the dipolar atom or molecule, and D is its 
dipole moment. The value of td can be compared with the average inter-particle spacing, 
r s ~ where p is the number density of the condensate and m the dimensionality. For 

rfj -C r s , the DDI is weak; in general, other contributions to the interaction, such as the 
s-wave scattering length a, will dominate (except if a is tuned to a sufficiently small value 2 ). 
For rr, > r s , the DDI will be the dominant interaction. The magnetic DDI is usually neg¬ 
ligible, only for the handful of atoms mentioned above, its effect has been observed, but it 
is difficult to increase the density such that rc > r s . Compared with the magnetic dipole 
moment of atoms, the electric dipole moment of heteronuclear molecules can be orders of 
magnitude larger, leading to large values for r^ (e.g. ru — 5 x 10 5 A for a fully polarized 
NaCs). Association of two atoms using a Feshbach resonance and transfer to the rovibra- 
tional ground state has been achieved for example for 'Li 133 Cs^, 40 K 8 'Rb^, 41 K 8 'Rb 10 and 
85 Rb 133 Cs n 12 . But it remains a challenge to produce a degenerate quantum gas of dipolar 
molecules. 

The anisotropy of the DDI leads to a measurable anisotropy of the speed of sound 1 ^, 
but also an anisotropic superfluid response® has been predicted. The attractive part of 
the DDI can give rise to roton or roton-like excitations in a dipolar Bose gas layei^ 15 19 . An 
anisotropic 2D quantum gas can be realized by tilting the polarization dipoles in a deep trap, 
and a stripe phase can form spontaneously 20 21 . For r£> r s , dipoles will crystallize without 
imposing an optical lattice 22 24 . A bilayered dipolar Bose gas can dimerize if the polarization 
direction in the two layers is the same 25 . Also the case of antiparallel polarization in two 
layers has been studied, where dipoles are perpendicular to the layer, but the orientation of 
the dipoles in one layer is opposite to that in the other layer®! 

In this work we study such a bilayer of bosonic dipoles with antiparallel polarization. 
The key result is that it is a self-bound liquid, unlike bilayers with parallel polarization, and 
therefore does not need external pressure in the form of a trap potential to stay together. 
We show that the liquid nature is a consequence of the attractive part of the inter-layer 
DDI, which leads to cohesion due to “dipole bridges” that effectively act as a glue to bind 
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all particles together. For our calculations we use a variational many-body theory, the 
hypernetted-chain Euler-Lagrange method (HNC-EL), which includes pair correlations. For 
comparison and validation, we use path integral Monte Carlo (PIMC) simulations. 

Methodology. A ID optical lattice slices a BEC into quasi-2D layers separated by a 
distance d. Since the dipole length m can easily exceed the typical d value of about 500nm, 
the DDI interaction between dipoles in different layers can lead to appreciable coupling. 
We consider here two translationally invariant layers A and R, approximate each layer as 
two-dimensional, and assume no tunneling occurs. With these simplifications we get two 
coupled 2D systems, i.e. a binary Bose mixture. The particles in the two layers shall be the 
same molecules, hence with same mass and dipole moment. In units of dipole length ro and 
the dipole energy Ed = h 2 /(mr 2 D ), the Hamiltonian is 



a and /3 index the layer, a,/3 G {A,B}, and i the particles within a layer. The primed sum 
indicates that for a = (3 we only sum over i ^ j, v a> p(\r it0l — r^j) is the DDI, in units of Ed, 
between dipole i at v i)0l in layer a and dipole j at r Jt rj in layer f3. We neglect short-ranged 
interactions compared to the DDI. The intralayer interaction (a = /3) is purely repulsive, 
v a ,a( r ) = 1/r 3 . The interlayer interaction, a ^ /3, is vab{ t ) = (2 d 2 — r 2 )/{d 2 -\ -r 2 ) 5 / 2 , which 
is repulsive for small r, but attractive for large r, and has a minimum at r m ; n = 2d. Since 
the average interlayer interaction vanishes, J d 2 rvAB{j" ) = 0, the coupling between layers in 
the ground state would vanish in a mean field approximation and the ground state energy 
would just be the sum of the energies of each layer. 

For the many-body ground state we use the variational Jastrow-Feenberg ansat^® con¬ 
sisting of a product of pair correlation functions for a multi-component Bose system, T 0 = 
ex P [\Y^a,pY!ij u a,p{Vi,a ~ r^l)]. Higher order correlations u a ^ n ( r ij0l ,r j}/3 ,r k}7 ) could be 
included, as is routinely done for single-component calculations. Past experience has shown 
that triplet correlations improve the ground state energy, leading to results very close to ex¬ 
act QMC simulations, but they do not change the qualitative picture. We therefore restrict 
ourselves to pair correlations, but check the results against PIMC simulations. The pair cor¬ 
relations u a fi{r) are determined from Ritz’ variational principle, i.e. from the coupled Euler- 
Lagrange equations, d(^o\H\^> 0 )/5u a} i3(r) = 0. They are solved by expressing u a ${r) in 
terms of the pair distribution function g a ^ (|r a —rg|) = j dr a ^ |T 0 | 2 where the inte- 
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gral is over all particles except one in layer a and one in layer f3, and p a = N a /V is the partial 
density of component a. u a A r ) and g a A r ) are related via the hypernetted-chain equations, 
and the equations to be solved for g a A r ) are the hypernetted-chain Euler-Lagrange (HNC- 
EL) equations. Certain “elementary” diagrams cannot be summed up exactly, and have to 
be approximated. We simply neglect them completely (HNC-EL/0 approximation); what 
we said about neglecting triplet correlations also applies to neglecting elementary diagrams. 
Details about the HNC-EL method can be found in Refs. 28 -^, and particularly for Bose 
mixtures in Refs.® 32 . The ffNC-EL/0 equations for an arbitrary number of components 
are (bold-faced capital letters denote matrices): 


W (k) 

V&(r) 


-- [S(fc)T(fc) + T(k)S(k) - 3T (k) 
S-\k)T(k)S-\k) 


g«A r ) v <*A r ) + 


h 2 

2 m a fi 



2 


+ (9aA r ) ~ ^Wa,p( T ) 


V ph (k ) = S _1 (fc)T(fc)S ~\k) - T (k) 


where m a ^ is the reduced mass (for the symmetric bilayer, h 2 /2m a = h 2 /2mp = 1/2 in 
dipole units). S a: p(k) is the static structure function S a: p(k) = 5 a p + y/p a pp FT [g a ,p ~ 1 ] 
(FT denotes Fourier transformation). The kinetic energy matrix, T a Ak) = S a p(h 2 k 2 /Am a: p), 
becomes T a Ak ) = <5 a /3w in our case. The ffNC-EL/0 equations can be solved iteratively. 
Usually the convergence is stable and fast, but close to an instability like the spinodal point 
discussed below, we use linear mixing between iterations to ensure convergence. 

Results. We calculated the ground state energy per particle, E(p)/N, as function of total 
density p = pa + Pb for different layer distances d. The interlayer DDf scales with d~ 3 , 
therefore the energy per particle, E/N , varies over a wide range, as can be seen in the top 
panel of Fig. [I] that shows E{p)/N for four values of d. A key result is that E(p)/N has a 
minimum at a certain equilibrium density p eq (d), where the pressure p vanishes: without an 
externally applied pressure provided e.g. by a trap potential, the total density of the bilayer 
system will adjust itself to p eq (d). Rather than expanding like a gas, a dipolar bilayer system 
with antiparallel polarization is a self-bound liquid. Despite the purely repulsive intralayer 
interaction, the partly attractive interlayer interaction provides the “glue” that binds the 
system to a liquid. The phenomenon of an effective intralayer attraction, mediated by 
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FIG. 1. (Color online) Top panel: Ground state energy per particle E/N versus the total density p , 
for several layer distances d. The results from HNC-EL/0 are shown as lines (small circles indicating 
the equilibrium density p eq ), the open symbols are from corresponding PIMC simulations. The 
filled symbol indicates E/N and p eq estimated from a PIMC simulation of a self-bound puddle of 
50 dipoles in each layer, see also 34 . Bottom panel: equilibrium density p eq versus d. 

particles in the other layer, is discussed in more detail below. 

During the iterative numerical optimization, convergence becomes very sensitive as the 
density p or the distance d between layers A and B is decreased, until the HNC-EL/0 
equations eventually fail to converge. Past experience with HNC-EL/0 is that a numer¬ 
ical instability usually has a physical reason. Indeed, as we will show below, there is a 
spinodal point where the homogeneous phase assumed in our calculation becomes unstable 
against phase separation by nucleation of puddles. Thus the equation of state E(p)/N for 
a homogeneous phase indeed ends at a critical density. 

Also shown in the top panel of Fig. [I] are the energies obtained with PIMC simulations. 
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The temperature is set to T = 0.5 (T = 0.25 for the smallest p), which is low enough 
that the thermal effect on E/N is smaller than the symbol size. The open triangles are 
bulk simulations of Na = Nb = 50 dipoles with periodic boundary conditions. The HNC- 
EL/0 results are upper bounds on E/N , consistent with a variational approach. The overall 
dependence of E/N on p and d is reproduced quite well with the HNC-EL/0 method, which 
is orders of magnitude faster than PIMC simulations. The black triangle shows the energy 
from a PIMC simulation of IV4 = Nb = 50 dipoles and layer separation d = 0.1 without 
periodic boundary conditions. Due to the liquid nature of the bilayer, the dipoles indeed 
coalesce into a puddle of finite density, given by p eq (d) apart from corrections due to the 
surface line tension. The density corresponding to the filled triangle is obtained from the 
radial density profile p(r) at r = 0 (set®) where r is defined relative to the center of 
mass of the puddle. Thermal evaporation was suppressed by choosing a low temperature of 
T = 1/16. Although this simulation of a finite cluster is not equivalent to bulk PIMC or 
HNC-EL/0 calculations, the central density of the puddle is close to p eq (<i) from HNC-EL/0. 

The equilibrium density p eq as function of d is shown in the bottom panel of Fig. 0 Peq{d) 
decreases rapidly with increasing d. For smaller d, the decrease is approximately p e q ~ d ~ 3 
and for larger d it is closer to p eq ~ d~ A . Based purely on the interlayer DDI Vab(t'), one 
would expect a scaling of p eq with the inverse square of r min = 2 d, leading to a scaling 
dr 2 . The deviation from d~ 2 is due to the kinetic energy. Only for very small d (very deep 
vab (r)), this simple picture can be expected to be valid, and indeed the p eq -curve becomes 
less steep for smaller d in the double logarithmic representation in Fig. [l] In this small-d 
regime of extremely strong interlayer correlation the HNC-EL method would not be reliable 
anymore. 

In Figj2] we show the intralayer and interlayer pair distributions, gAA{f) and Qab/i"), 
in the top and middle panel for progressively smaller layer distance d up to the smallest 
numerically stable value d = 0.063, for p — 1. The growth of a strong correlation peak 
in (]ar(j) as d is decreased is a direct consequence of the increasingly deep attractive well 
of vab(j" ) around r min = 2d. But also gAA( r ) develops additional correlations, seen as a 
shoulder in the top panel. The additional correlations are best seen in the difference to 
the uncoupled (d = 00) limit 5X4 ( r )> ^Qaa( r ) = gAA(j') — 5X4 ( r )> shown in the bottom 
panel. This additional positive correlation between dipoles in the same layer is mediated by 
dipoles in the other layer: the attraction between a dipole in layer A and a dipole in layer 
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FIG. 2. (Color online) Top panel: Intralayer pair distributions gAAir) at density pa = Pb = 0.5 for 
progressively smaller layer distance d. Middle panel: Corresponding interlayer pair distributions 
9 AB{r), with circle indicating the maxima of gAB{f)- Bottom panel: incremental intralayer pair 
distributions AgAA ^’) = gAA{f) — 5 ) 4 ( 4 ( r )> he. change from uncoupled layers. The inset in the 
middle panel shows the positions of the maxima of gAA^r) and 5,4 n(r), respectively, as function of 
d. The inset in the bottom panel sketches the attractive forces between dipoles in different layers. 

B, that leads to a peak at distance r m - a bit larger than 2d due to zero-point motion 
induces an effective attraction between the dipole in A and another dipole in A, leading to 
peak at about twice the distance, 2r m . The inset in the middle panel shows the maxima of 
9 ab(t) and A gAA{f) (indicated by circles in the plots of gAB^r) and A qaa{t)) as function 
of distance d. Indeed the peaks of AgAA(r ) are located at about twice the distance of the 
peaks of gAB{r )• This effective intralayer attraction, induced by the interlayer attraction, is 
illustrated by a simple ID picture in the inset in the bottom panel, which also illustrates 
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a preference for a certain interparticle spacing, i.e. density, where “dipole bridges” (blue 
lines) can form. The present 2D situation is more complicated, but our results for the pair 
correlations demonstrate this picture is approximately valid. 



FIG. 3. (Color online) Speed of density fluctuations c\ and speed of concentration fluctuations C 2 , 
as function of layer distance d for three densities, p = 0.1; 1; 10 (top, middle, and bottom panels). 
Full lines show the Bijl-Feynman approximation and symbols are the thermodynamic estimates. 

The identification of the low density instability with a spinodal point can be proven by 
calculating the long wavelength modes. For two coupled layers, there are two excitation 
modes, e\ o for any given wave number k, a density mode and a concentration mode. In the 
long wavelength limit, k —>• 0, each mode can be characterized by the speed of a density or 
concentration fluctuation, C\ and C 2 , respectively. At the spinodal point, C\ vanishes which 
means that the system becomes unstable against infinitesimal k —s- 0 fluctuations of the 
total density, triggering the spinodal decomposition. The easiest way to calculate c\ and 
C 2 is the Bijl-Feynman approximation (BFA) for the excitation energies 6j(/c), i.e. solving 
the generalized eigenvalue problem ^0 = eS(fc)0. For strong correlations, the BFA gives 
only a rough idea of the true excitation structure, e.g. the BFA for the roton energy of 












superfluid 4 He is off by a factor of two. However, it describes the low momentum limit of 
the dispersion relation very well, which is what we need for q. For a symmetric bilayer, 
the eigenvalues are t\ j 2 (fc) = -^-( Saa {k) ± S A B(k ))~ 4 and the associated eigenvectors are 
0i ~ (1,1) and 0 2 ~ (1,-1). 0i describes total density fluctuations where particles in 
different layers move in phase and 02 describes concentration fluctuations, where particles 
in different layers move out of phase, and the total density is constant. For small k , ei (k) < 
62 (k), i.e. the density mode has lower energy than the concentration mode. For k 0 we get 
c i,2 = \{S'aa i S'ab) \ w h ere we abbreviated the derivatives S' a * = dS a ^(k)/dk\k=o- For 
single-component Bose systems, it is known that long-wave length limit of S(k) obtained 
with HNC-EL are biased by the approximation made for elementary diagrams (omitted 
here altogether). This leads to an inconsistency between the speed of sound c obtained from 
the HNC-EL approximation for S(k) and the thermodynamic relation between c and the 
energy, c 2 = (in dipole units). In order to assess the reliability of our results for 

Ci and c 2 , we compare the BFA values obtained from S a ^, Ci >2 = \{S AA ± S' AB )~ l , with the 
generalization of the thermodynamic relation between q i2 and the energy to binary systems, 
cl 2 — p( e AA T 6 A b) / 2 where e a p is the second derivative of E/N with respect to p a and 

pP 

In Fig. [3] the results for Ci and c 2 obtained with the two methods are shown as function 
of layer distance d for densities p = 0.1; 1; 10. As the coupling between layers is increased 
by reducing d, ci and c 2 behave differently. The speed of concentration fluctuations c 2 
increase (without actually diverging) while the quantity of main interest, the speed of density 
fluctuations ci decreases to zero, in agreement with the interpretation of the instability 
as a spinodal point. The critical distance where Ci vanishes is lower for higher p, hence 
increasing the density for a given d makes the system more stable. The BFA for ci )2 and 
their thermodynamic estimates agree qualitatively, but differ especially for the interesting 
regime near the spinodal point where Ci —y 0. The Bijl-Feynman values for ci appear to go 
to zero linearly and at slightly smaller d, while the thermodynamic estimates approach zero 
more steeply, possible in a non-analytic fashion. Unfortunately, these uncertainties preclude 
a meaningful analysis of critical exponents for Ci(d) or ci(p). Monte Carlo simulations, 
including a finite size scaling analysis, may shed more light on this question, but would 
certainly require very large simulations, beyond the scope of this paper. 

In conclusion, we have shown that a dipolar bilayer with antiparallel polarization in 
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the two layers constitute a self-bound liquid, evidenced by a minimum of E(p)/N at a 
finite density. The bilayer relaxes to a stable equilibrium at a finite total density and 
requires no external pressure coming from a trapping potential, which makes is possible 
to study homogeneous quantum phases. Comparison with exact PIMC simulations shows 
good agreement with results obtained with the HNC-EL/0 method. As expected for a liquid, 
the equation of state ends at a spinodal point where the bilayer becomes unstable against 
infinitesimal long-wavelength perturbations, thus the speed of total density fluctuations 
approaches zero. Finite size PIMC simulations confirm that the system indeed coalesces 
into a puddle with a flat density profile given by the equilibrium density 34 . 
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Supplement 


Comparison with PIMC 

An important question is how accurate are our ground state results obtained with the 
variational HNC-EL/0 method, where we neglect elementary diagrams and higher than two- 
body correlations in the ansatz for the wave function. We performed path integral Monte 
Carlo (PIMC) simulations to assess the qualtity of our HNC-EL/0 results and found good 
agreement between HNC-EL/0 energies and PIMC energies at low temperature. In this 
supplement we present additional comparisons of the static structure matrix, as well as 
results for finite systems, where self-bound “puddles” are formed. All simulations were done 
with 50 particles per layer; bulk simulations used quadratic simulation boxes with periodic 
boundary conditions with a box size adjusted to achieve a given density. For antiparallel 
bilayers, cut-off corrections to the dipole-dipole interaction cancel each other and therefore 
are not needed. 

For predicting the spinodal instability, but also for calculating excitation properties, an 
important quantity is the static structure matrix, S a p(k). In Figj4]we compare Saa( k) and 
Sab( k) obtained with HNC-EL/0 (lines) and PIMC (symbols), at a total density of p = 1. 
The temperature in the PIMC simulation was T = 0.5 which was low enough that S a p{k) 
did not change upon lowering the temperature. We see that the HNC-EL/0 approximation 
works well, considering that intralayer correlations are quite strong. For d = 0.3 and d = 0.15 
PIMC simulations predicts slightly more pronounced peaks and troughs, but HNC-EL/0 
calculations are faster by several orders of magnitude. For d = 0.1, the agreement is also 
very good, except for the two smallest k values possible in a simulation box of side length 10, 
k = 27t/ 10 ~ 0.63. In the PIMC results, both Saa{}S) and Sab{}S) turn up sharply for this 
smallest k value. When we reduce d even more, this apparent peak at k — 0 grows very large. 
This peak has a very simple reason: as we approach the spinodal point by reducing d, we 
enter the metastable regime of the phase diagram, where E/N as function of total density 
has a negative slope. In this regime a finite perturbation can lead to a collapse and the 
system phase separates. Since there is no trial wave function in PIMC that could prevent 
that, this collaps indeed happens as we go to far below the equilibrium density. Indeed, 
Monte Carlo snapshots such as in Fig. [5] show density fluctuations already for d = 0.1 that 
resemble small “bubbles”. In Fig. [5] red and blue dots, connected by lines, are the beads of 
the discretized imaginary time paths sampled in PIMC; each bead is a particle at a discrete 
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FIG. 4. (Color online) Static structure functions SAA{k ) (upper panel) and SAB(k) (lower panel) 
obtained with HNC-EL/0 (lines) for a total density p = 1 and progressively smaller layer distances 
d as indicated in the upper panel. The symbols show SAA(k) and SAb(&) obtained by PIMC 
simulations at T = 0.5K for distances down to d = 0.1. 

time step. For even lower d or lower total densities we observe a clear decomposition into a 
droplet and a low-density gas, see next paragraph. A large peak at k = 0 can then be seen 
as a zero momentum Bragg peak due to phase separation. 

As final confirmation of the liquid nature of dipolar bilayers with antiparallel polarization, 
we show the results of a PIMC simulation without periodic boundary conditions and without 
any external trap potential in the planar direction. A two-dimensional gas would of course 
spread out indefinitely, while a liquid will coalesce into a droplet of finite density. For 
a large enough droplet such that effects of surface line tension are negligible, the density 
inside the droplet is given by the equilibrium density p eq , i.e. the density of a bulk system 
at zero pressure. In Fig|6]we show the radial density profile p(r ) for 50 dipoles in each layer 
separated by d = 0.1, where r is measured relative to the center of mass. In order to prevent 
evaporation we set the temperature to T = A. p( r ) i s approximately constant for r < 4 and 
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FIG. 5. (Color online) PIMC simulations snapshot for p = 1 and d = 0.1 (T = 0.5). Red and 
blue chains are the dipoles in layer A and B. 

quickly falls to zero for larger r. This is the behavior expected for the density profile of a 
self-bound liquid, and very different from the density profile of a quantum gas in a trap. 
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FIG. 6. (Color online) Density profile p(r) for a self-bound droplet of 50 dipoles in each layer. 
Inside the droplet the density is approximately constant, with a value close to the equilibrium 
density p eq of a bulk system at zero pressure. The distance is d = 0.1 and the temperature was 
set to T = j|r, which is low enough to prevent evaporation. The inset shows a snapshot of the 
simulation, with red and blue indicating the dipoles in the layers A and B, respectively, at the 
imaginary time steps of the paths of PIMC. 
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